{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nTkZ0i6Yzhym",
    "papermill": {
     "duration": 0.007255,
     "end_time": "2021-01-04T05:10:58.373630",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.366375",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Vaccine RCT Examples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5MneoLFoegqr"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wGGzchhjzhyy",
    "papermill": {
     "duration": 0.006053,
     "end_time": "2021-01-04T05:10:58.386521",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.380468",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "  \n",
    "# Polio RCT\n",
    "\n",
    "One of the earliest randomized experiments were the Polio vaccination trials conducted by the Public Health Service in 1954.  The question was whether Salk vaccine prevented polio.  Children in the study were randomly assigned either a treatment (polio vaccine shot) or a placebo (saline solution shot), without knowing which one they received. The doctors in the study, making the diagnosis, did not know whether a child received a vaccine or not. In other words, the trial was a double-blind, randomized control trial.  The trial had to be large, because the rate at which Polio occured in the population was 50 per 100,000.  The treatment group saw 33 polio cases per 200,745; the control group saw 115 cases per 201,229. The estimated average treatment effect is about\n",
    "$$\n",
    "-40\n",
    "$$\n",
    "with 95% confidence band (based on approximate normality of the two sample means and their differences)\n",
    "$$[-52, -28].$$\n",
    "The confidence suggests that the Polio vaccine **caused** the reduction in the risk of polio.\n",
    "\n",
    "The interesting thing here is that we don't need the underlying individual data to evaluate the effectivess of the vaccine. This is because the outcomes are Bernoulli random variabales, and we have enough information to compute the estimate of ATE as well as the confidence intervals.\n",
    "\n",
    "\n",
    "We also compute Vaccine Efficacy metric, which (I googled  [CDC](https://www.cdc.gov/csels/dsepd/ss1978/lesson3/section6.html)) refers to the following measure:\n",
    "$$\n",
    "\\operatorname{VE} = \\frac{\\text{Risk for Unvaccianted - Risk for Vaccinated}}{\\text{Risk for Unvaccianted}}.\n",
    "$$\n",
    "It describes the relative reduction in risk caused by vaccination.\n",
    "\n",
    "\n",
    "It is staighborward to get VE estimate by just plugging-in the numbers, but how do we get the approximate variance estimate. I am too lazy to do calculations for the delta methods, so I will just use a simulation (a form of approximate bootstrap) to obtain the confidence intervals.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "id": "n0AUYEIczhy2",
    "papermill": {
     "duration": 1.94521,
     "end_time": "2021-01-04T05:11:00.339761",
     "exception": false,
     "start_time": "2021-01-04T05:10:58.394551",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Here we calculate the overall effectiveness of the vaccine and construct confidence intervals for it\n",
    "NV = 200745  # number of vaccinated (treated)\n",
    "NU = 201229  # number of unvaccinated (control)\n",
    "RV = 33 / NV  # average outcome for vaccinated\n",
    "RU = 115 / NU  # average outcome for unvaccinated\n",
    "VE = (RU - RV) / RU  # vaccine efficacy\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "\n",
    "print(f\"Incidence per 100000 among treated: {RV * 100000:.4f}\")\n",
    "print(f\"Incidence per 100000 among controlled: {RU * 100000:.4f}\")\n",
    "\n",
    "# treatment effect: estimated reduction in incidence per 100k people\n",
    "effect = 100000 * (RV - RU)\n",
    "\n",
    "print(f\"Estimate of effect on occurances per 100,000: {effect:.4f}\")\n",
    "\n",
    "# variance, standard deviation and confidence interval of ATE using\n",
    "# the fact that outcomes are Bernoulli\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "var_effect = (100000**2) * (var_rv + var_ru)\n",
    "std_effect = np.sqrt(var_effect)\n",
    "\n",
    "print(f\"Standard error for ATE is: {std_effect}\")\n",
    "\n",
    "ci_effect = [effect - 1.96 * std_effect, effect + 1.96 * std_effect]\n",
    "\n",
    "print(f\"95% confidence interval is [{ci_effect[0]:.4f}, {ci_effect[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nY4oaYiIegq5"
   },
   "outputs": [],
   "source": [
    "# we use approximate bootstrap to calculate confidence\n",
    "# intervals for vaccine efficacy via Monte Carlo draws\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV + np.random.normal(0, 1, B) * np.sqrt(var_rv)\n",
    "RUs = RU + np.random.normal(0, 1, B) * np.sqrt(var_ru)\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "11nKcV5pegq7"
   },
   "outputs": [],
   "source": [
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "79c7e3d0-c299-4dcb-8224-4455121ee9b0",
    "_uuid": "d629ff2d2480ee46fbb7e2d37f6b5fab8052498a",
    "id": "In0Vkgi1zhy4",
    "papermill": {
     "duration": 0.007763,
     "end_time": "2021-01-04T05:11:00.355927",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.348164",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Pfizer/BNTX Covid-19 RCT\n",
    "\n",
    "Here is a link to the FDA [briefing](https://www.fda.gov/media/144245/download) and an interesting [discussion](\n",
    "https://garycornell.com/2020/12/09/statistics-in-the-pfizer-data-how-good-is-the-vaccine/?fbclid=IwAR282lS0Vl3tWmicQDDhIJAQCMO8NIsCXyWbUWwTtPuKcnuJ2v0VWXRDQac), as well as data.\n",
    "\n",
    "Pfizer/BNTX is the first vaccine approved for emergency use to reduce the risk of Covid-19. Volunteers were randomly assigned to receive either a treatment (2-dose vaccination) or a placebo, without knowing which they received. The doctors making the diagnoses did not know now whether a given volunteer received a vaccination or not. The results of the study are given in the following table:\n",
    "\n",
    "![pfizer/biontech result table](https://lh6.googleusercontent.com/oiO6gYom1UZyrOhgpFx2iq8ike979u3805JHiVygP-Efh1Yaz2ttyPcgWKlT1AqHDM4v46th3EPIkOvRLyXA0fNUloPL-mL9eOFmSAzfbNOHyCZSQ0DyzMhcFUtQuZ520R5Qd2lj)\n",
    "\n",
    "Here we see both the overall effects and the effects by age group. The confidence intervals for the overall ATE are tight and suggest high effectivness of the vaccine. The confidence intervals for the age groups 65-74 and 75+ are much wider.  We could group the 65-74 and >75 groups to evaluate the effectiveness of the vaccine and also narrow down the width of the confidence band.\n",
    "\n",
    "In this case, the reported results are for vaccine effectiveness. We use the same approach as above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wLulUAGdzhy5",
    "papermill": {
     "duration": 0.007593,
     "end_time": "2021-01-04T05:11:00.371394",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.363801",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the code cell below  we calculate the overall effectiveness of the vaccine and construct confidence intervals for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fmNIIhxuzhy9",
    "papermill": {
     "duration": 0.497816,
     "end_time": "2021-01-04T05:11:00.876971",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.379155",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "NV = 19965\n",
    "NU = 20172\n",
    "RV = 9 / NV\n",
    "RU = 169 / NU\n",
    "VE = (RU - RV) / RU\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "\n",
    "print(f\"Incidence per 100000 among treated: {RV * 100000:.4f}\")\n",
    "print(f\"Incidence per 100000 among controlled: {RU * 100000:.4f}\")\n",
    "\n",
    "# treatment effect: estimated reduction in incidence per 100k people\n",
    "effect = 100000 * (RV - RU)\n",
    "\n",
    "print(f\"Estimate of effect on occurances per 100,000: {effect:.4f}\")\n",
    "\n",
    "# variance, standard deviation and confidence interval of ATE using\n",
    "# the fact that outcomes are Bernoulli\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "var_effect = (100000**2) * (var_rv + var_ru)\n",
    "std_effect = np.sqrt(var_effect)\n",
    "\n",
    "print(f\"Standard error for ATE is: {std_effect}\")\n",
    "\n",
    "ci_effect = [effect - 1.96 * std_effect, effect + 1.96 * std_effect]\n",
    "\n",
    "print(f\"95% confidence interval is [{ci_effect[0]:.4f}, {ci_effect[1]:.4f}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "3coxMaz5egrG"
   },
   "outputs": [],
   "source": [
    "# we use approximate bootstrap to calculate confidence\n",
    "# intervals for vaccine efficacy via Monte Carlo draws\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV + np.random.normal(0, 1, B) * np.sqrt(var_rv)\n",
    "RUs = RU + np.random.normal(0, 1, B) * np.sqrt(var_ru)\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sgE4OHh5zhy-",
    "papermill": {
     "duration": 0.009208,
     "end_time": "2021-01-04T05:11:00.895842",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.886634",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the code cell below, we calculate the effectiveness of the vaccine for the two groups that are 65 or older"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pvzbZxlAzhy_",
    "papermill": {
     "duration": 0.477817,
     "end_time": "2021-01-04T05:11:01.383093",
     "exception": false,
     "start_time": "2021-01-04T05:11:00.905276",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "NV = 3239 + 805\n",
    "NU = 3255 + 812\n",
    "RV = 1 / NV\n",
    "RU = (14 + 5) / NU\n",
    "VE = (RU - RV) / RU\n",
    "\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = RV + np.random.normal(0, 1, B) * np.sqrt(var_rv) + 1e-10\n",
    "RUs = RU + np.random.normal(0, 1, B) * np.sqrt(var_ru) + 1e-10\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"Two-sided 95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.05))\n",
    "\n",
    "print(f\"One-sided 95% confidence interval of VE is [{CI_VE:.4f}, 1]\")\n",
    "\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UQw0zHWGegrK"
   },
   "source": [
    "Let's try the Binomial parametric bootstrap using the fact that outcomes are Bernoulli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "W-nfdQ9UzhzA",
    "papermill": {
     "duration": 0.539729,
     "end_time": "2021-01-04T05:11:01.934674",
     "exception": false,
     "start_time": "2021-01-04T05:11:01.394945",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "NV = 3239 + 805\n",
    "NU = 3255 + 812\n",
    "RV = 1 / NV\n",
    "RU = (14 + 5) / NU\n",
    "VE = (RU - RV) / RU\n",
    "\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {VE:.4f}\")\n",
    "var_rv = RV * (1 - RV) / NV\n",
    "var_ru = RU * (1 - RU) / NU\n",
    "\n",
    "np.random.seed(123)\n",
    "B = 10000\n",
    "RVs = np.random.binomial(NV, RV, size=B) + 1e-10\n",
    "RUs = np.random.binomial(NU, RU, size=B) + 1e-10\n",
    "VEs = (RUs - RVs) / RUs\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.025, .975))\n",
    "\n",
    "print(f\"Two-sided 95% confidence interval of VE is [{CI_VE[0]:.4f}, {CI_VE[1]:.4f}]\")\n",
    "\n",
    "CI_VE = np.quantile(VEs, (.05))\n",
    "\n",
    "print(f\"One-sided 95% confidence interval of VE is [{CI_VE:.4f}, 1]\")\n",
    "\n",
    "\n",
    "sns.kdeplot(VEs, shade=True)\n",
    "plt.title(\"Approximate distribution of VE estimates\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jPZWoSzsegrM"
   },
   "source": [
    "# More Accurate Confidence Intervals\n",
    "\n",
    "The table from the FDA uses exact inference by inverting tests based on the exact binomial nature of the outcome variable. This method is known as the Cornfield Procedure to find the exact confidence interval on the estimate of vaccine efficacy.\n",
    "\n",
    "Another typical more accurate approximation is to approximate the log of the risks by a normal distribution, construct confidence intervals for the quantity $\\log(RV) - \\log(RU)$ and then invert. This is for instance method C in [this paper](https://www.jstor.org/stable/pdf/2530610.pdf) and is implemented in the statsmodels package (and also in the `scipy.stats.contengency_table` package)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NaMX-TVPegrN"
   },
   "outputs": [],
   "source": [
    "tb = sm.stats.Table2x2([[9, 19965 - 9], [169, 20172 - 169]])\n",
    "tb.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d-w2izjvegrO"
   },
   "outputs": [],
   "source": [
    "print(f\"Overall vaccine efficacy is {1 - tb.riskratio:.4f}\")\n",
    "lb, ub = tb.riskratio_confint(alpha=.05)\n",
    "print(f\"Approximate 95% confidence interval of VE: [{1 - ub:.4f}, {1 - lb:.4f}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d9rPuODEegrP"
   },
   "source": [
    "We see that this matches much more accurately the table from the FDA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8ntswBDwrfTt"
   },
   "source": [
    "Next we repeat this procedure to find a more accurate confidence interval on vaccine effectiveness for the two groups that are 65 or older. Here we see a big discrepancy between various asymptotic approaches and the exact finite-sample inference. This occurs because the binomial counts are too low for central limit theorems to work successfully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8GaXZidxrgsr"
   },
   "outputs": [],
   "source": [
    "# Exact CI exploiting Bernoulli outcome for the two groups that are 65 or older\n",
    "NV = 3239 + 805\n",
    "NU = 3255 + 812\n",
    "RV = 1 / NV\n",
    "RU = (14 + 5) / NU\n",
    "VE = (RU - RV) / RU\n",
    "\n",
    "tb = sm.stats.Table2x2([[1, NV - 1], [19, NU - 19]])\n",
    "tb.summary()\n",
    "\n",
    "print(f\"Overall vaccine efficacy is {1 - tb.riskratio:.4f}\")\n",
    "lb, ub = tb.riskratio_confint(alpha=.05)\n",
    "print(f\"Approximate 95% confidence interval of VE: [{1 - ub:.4f}, {1 - lb:.4f}]\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "papermill": {
   "duration": 8.643307,
   "end_time": "2021-01-04T05:11:02.054671",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-01-04T05:10:53.411364",
   "version": "2.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
